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Abstract 

Domain-wall fermions preserve chiral symmetry up to terms that de- 
crease exponentially when the lattice size in the fifth dimension is taken 
to infinity. The associated rates of convergence are given by the low-lying 
eigenvalues of a simple local operator in four dimensions. These can be 
computed using the Ritz functional technique and it turns out that the 
convergence tends to be extremely slow in the range of lattice spacings rel- 
evant to large-volume numerical simulations of lattice QCD. Two methods 
to improve on this situation are discussed. 
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Introduction 



The idea to realize 4- dimensional (4D) chiral fermions on the lattice by coupling 
5D fermions to a 4D domain wall Jl| has attracted a lot of attention in the lattice 
community (for a review see [§]). Although originally designed to construct chiral 
gauge theories, the idea can also be applied to lattice QCD in order to preserve 
the global chiral symmetry at zero quark mass 0, ||. In this 5D Wilson- 

Dirac operator is chosen with N slices in the extra dimension and appropriate 
Dirichlet boundary conditions in the fifth dimension. In the limit N — > oo, chiral 
zero modes exist as surface modes on the 4D boundaries, even at finite lattice 
spacing. 

The bulk fermionic degrees of freedom are massive and can be shown to decouple 
in the continuum limit: the action of the 5D system is equivalent to the one 
corresponding to a 4D Dirac operator describing the boundary chiral modes H; 
similarly, the propagator of the boundary fields can be obtained from the inverse 
of the same 4D operator ||. The chirality of the surface modes in the limit 
iV — > oo then follows || from the fact that, in this limit, this 4D Dirac operator 
satisfies the Ginsparg- Wilson (GW) relation 0-0, which implies an exact lattice 
chiral symmetry ||10|| . Thus the 5D domain- wall construction in the limit iV — > oo 
is completely equivalent to a 4D lattice formulation of Ginsparg- Wilson fermions, 



and satisfies all the properties that follow from the exact chiral symmetry [pi 11 



12, 10, 13 . Moreover, if the continuum limit is taken in the extra dimension, this 



4D formulation coincides with that using Neuberger's fermion operator [14 



The introduction of an extra dimension makes domain-wall fermions more de- 
manding numerically than standard Wilson fermions (the equivalent 4D formu- 
lation is similarly more demanding, owing to the non-ultralocality of the action). 
Nevertheless the advantage of preserving an exact chiral symmetry might com- 
pensate for the higher cost in some cases. An analysis in the free theory showed 
that the convergence to the exact operator as a function of N is rapid 0, |I5| . 
This gave rise to the hope that also in the interacting case domain-wall fermions 
could be used without too much computational overhead. However, in realistic 
simulations, there are indications that the convergence rate deteriorates rapidly 
at large values of the gauge coupling, and much larger values of iV are indeed 



needed |T6|-[jT8|, leading to a substantial computational cost. 
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In a recent paper [19|| , the problems found in practical simulations were traced 
back to the appearance of very small eigenvalues of a certain 4D operator, which 
control the rate of convergence in N. We have performed an independent study 



and confirm the analysis in flTJ| . In addition, we discuss a new method to improve 
the domain-wall fermion operator, which differs from the one proposed in [|TIJ and 
proves to work better numerically. 



Five dimensional theory 



In this section we establish our notation and collect some useful formulae, the 
derivation of which can be found in P, ^ |5|, |6|, pOp. The 5D domain wall operator 
we consider here is defined as 

V = \ {75(5; + d.) - a s d* s d s } + M, (1) 

where s denotes a lattice site in the fifth direction (a s < s < a s N), a s is the 
corresponding lattice spacing, and d* and d s are the usual forward and backward 
derivatives. The operator M in eq. (|l|) is obtained from the standard 4D Wilson 
operator by 

M = D w -m (2) 

with 

^w = ^{7M(v; + v M )-av;v M }. (3) 

Here V* and V M are the gauge covariant forward and backward derivatives and 
a is the lattice spacing in the four physical dimensions /i = 1,...,4. The mass 
parameter mo obeys 

m > 0, a s mo < 2, am < 2. (4) 

Note that the lattice spacings a s and a can be different. The boundary conditions 
are fixed through 

P+V(0,x) = P.tfj(a s N + a s ,x) = , (5) 
where P± = |(1 ± 75). 
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Appropriate interpolating fields for the quarks constructed out of the left and 
right boundary modes are 



q(x) = P-ij){a a ,x) + P + 4j{Na s ,x) , (6) 
q(x) = $(a s , x)P + + ip(Na s , x)P-. (7) 

A mass term can then be introduced by adding to eq. ([!]) the term 

^m{tp(a s ,x)P + ifj(Na s ,x) +ifj(Na s ,x)P_ifj(a s ,x)} . (8) 

The two-point function of the quark fields is related to an effective 4D operator 
D N by [| 

(Q(x)q(x)) = 2 -^£E, (9) 

with 

D m ,N = (1 - 7^am)D N + m. (10) 
In terms of the operators K±, 

1 1 a s M . . 



D N is given by 



aD N = ^+^ K + N + RN - (12) 



From eq. (|T2|) it is straightforward to show that 



aD = lim aD N = 1 + -f 5 e(K + - KJ) , (13) 

N— >oo 



which can also be written as 20 



aD = 1-AiA^A)- 1 ' 2 (14) 
A = -a s M(2 + a s M)~ 1 . (15) 

It follows easily from this expression that D satisfies the Ginsparg-Wilson rela- 
tion, the only difference to Neuberger's operator being the different definition of 
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A. Indeed, Neuberger's operator is readily obtained from eqs. ([L4]) and (|T^) by 
taking the limit a s — > 0. 

Similarly, in the limit A" — > oo, the fermion determinant of the 5D formulation 
can be written in terms of the determinant of -D m ,jv, up to local subtractions. 
The 5D formulation is thus completely equivalent to a 4D lattice formulation of 
Ginsparg- Wilson fermions satisfying an exact chiral symmetry. 
A final necessary condition for this formulation to be an acceptable regularization 
of QCD is that the operator of eq. (Ill) is local. Indeed, it has been shown by 



Kikukawa that both the operators in eqs. ([H]) and (|T2|) are exponentially 
localized for smooth enough gauge fields, satisfying a plaquette bound |22|, |23[ . 
In realistic simulations of domain-wall fermions, A" is finite. In this situation, 
the chiral symmetry is broken by the residual terms SD = — D. It may be 
speculated that one could include a small additive quark mass renormalization, in 
order to get rid of these chirality breaking effects. This is, however, only justified 
by universality arguments if the subleading corrections in A" in the action are 



local. The result of |21| shows that this is indeed the case since SD is also local. 
However it is important to stress that the exponential localization of SD only 
sets in at distances of O(N). This can be shown already in the free case. On the 
other hand, in practical simulations the typical lattice sizes used are not much 
larger than A" and consequently SD is not local at the distances probed. In 
this situation, a quark mass renormalization cannot cancel the chirality-breaking 
effects induced by SD. 

The convergence rate in N 

For gauge field configurations with a restricted plaquette value, the operator A^A 



has a spectral gap [22 



< u < A ] A < v, (16) 

ensuring the exponential convergence in N of D^. The minimum rate of conver- 
gence is given by 

oo = mini^], uJi = In | + J^. , (17) 
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where Aj are the eigenvalues of A^A. 

However, in realistic simulations the plaquette bound is not satisfied and it is 
important to study the convergence rate u for the values of (3 and mo at which 
large-scale numerical simulations are performed nowadays. 

The eigenvalues A, which determine u, can be obtained through the generalized 
4D eigenvalue equation 

a 2 s M ] M^) = \(2 + a s M)\2 + a s M)^. (18) 

It is either the minimum or maximum eigenvalue of A' A that minimizes cu,. 
These eigenvalues can be obtained by minimizing (maximizing) the generalized 
Ritz functional 

(ip\a 2 s M^M\ilj) 



(V|(2 + a a M)t(2 + a.M)|V) 
using a straightforward generalization of the algorithm described in [24] []. 



(19) 



Eigenvalues above the lowest one can be computed by modifying the operator 
M^M in the numerator of eq. (|19|) in such a way that the already computed 
eigenvalues are shifted to larger values. For example, this can be achieved by 
substituting M^M by M^M + £\ .M^M^) (^[(2 + a s M)t(2 + a s M) with act = 
(1 — Aj)/Aj and Aj, ipi being the already computed eigenvalues and vectors. Notice 
that in this method no inversion of the matrix (2 + a s M)'(2 + a s M) is needed. 
We have studied numerically the convergence rate in the quenched approximation. 
We find that it is always controlled by the lowest eigenvalue A m i n of A^A. In Figs. |I] 
and we show the inverse convergence rates oj~ x corresponding to the five lowest 
eigenvalues of A^A at (5 = 5.85 on an 8 3 • 16 lattice, and at (3 — 6.0 on a 16 3 • 32 
lattice, respectively. In both cases we have set a s m Q = am = 1.8, which is a 
typical value used in previous simulations. 

Figures and |2| give a rather pessimistic view of the convergence of domain- 
wall fermions to the exact operator: they imply that several hundred or even 
thousand slices in the extra dimension would be needed to achieve a reasonable 
approximation. Clearly, this would render domain-wall fermions impracticable. 



Our results are consistent with the findings in [19 



1 A more detailed description of the algorithm can be obtained from the authors on request. 
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Figure 1: Inverse convergence rate to^ 1 for the five lowest eigenvalues of A^A 
(open symbols correspond to the lowest eigenvalue) as a function of Monte Carlo 
time £mc; at /3 = 5.85 on an 8 3 • 16 lattice. 
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;ure 2: The same as Fig. |T] at (3 = 6 on a 16 3 ■ 32 lattice. 
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It should be noted, however, that a very similar situation was found for Neu- 
berger's operator |f22| . Also in this case very small eigenvalues of the corre- 
sponding A* A = M^M occur, turning the computation of its inverse square root 
extremely costly. 

Acceleration of convergence 



In the case of Neuberger's operator, the bad convergence behaviour resulting 
from the low-lying eigenvalues of M^M could be cured by treating these modes 
exactly [^5, ^6|. It is natural to look for a similar trick also for domain- wall 



fermions, given the similarity of the two constructions. We have found two ways 



of achieving this. The first method is equivalent to the one described in |T9 
so we will not give any details here. The corresponding improved 5D operator 
differs from the standard one by boundary terms. 

We have tested the inversion of this improved operator, V impr , by solving the 
linear equation Di mpr X = Y for a given source Y. As a numerical solver we have 
used a conjugate gradient method with a stopping criterion ||r||/||X|| < e, where 
r = Vi mpr X — Y is the residual vector and e was set to e = 10~ 8 . We found 
that when using such a relatively low value of e the conjugate gradient method 
behaves very poorly: for a number of configurations at /3 — 5.85 on an 8 3 ■ 16 
lattice, the norm of the residual vector developed a very long tail at rather small 
values of ||r|| < O(10~ 5 ). This resulted in a very large number of iterations in 
the conjugate gradient algorithm before it converged to the desired accuracy. We 
suspect that some subtle cancellations occur in the improved operator leading to 
unexpectedly large rounding error effects. 

Since this behaviour of the conjugate gradient algorithm was rather unsatisfac- 
tory, we developed an alternative improvement method. The key observation for 



the new improved 5D operator is that the relations in eq. (15) and eq. (TH|) hold 
true for any choice of M as long as 

M f = 7 5 M 75 , det(2 + a s M) ^ . (20) 

This fact may be used to construct an improved M for which the very low eigen- 
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values of A^A disappear. A possible form of M that achieves this is given by 

r 

a s M = a s M - X ki™k ® w}j s , (21) 

k,l=l 

where 

^Av k = a k v k , k = 1, • • • , r, u { ) = 5 H . (22) 

and 

w fc = (2 + a s M)^v k . (23) 

Finally 

(X- 1 )^ = 25 kl (a k - ak)- 1 + (v fc , (2 + a s M) T5 ^) . (24) 

The corresponding 5D operator T> is given by eq. (|l|) after substituting M by M. 
Notice that the improved operator differs from the original one also in the bulk 
and not just at the boundary. 
After some algebra it can be shown that 

r 

A = -a s M(2 + a s Niy l = A + ^(d fc - a k )^v k <g> v\ . (25) 

k=l 

It is now easy to see that 75A has the same eigenvectors as 75A, but all eigenval- 
ues a k , k — 1, • • • ,r, are replaced by a k . The limit iV — > 00 of the correspond- 
ing improved 4D operator Dn is the same as that of the original Dn provided 
sign(dfc) = sign(afc). However, the approach to this limit is faster for D N if 
the lowest eigenvalues of A^A, X k = acj,, are larger than those of A^A, i.e. if 
\a k \ > \a k \. 

The concrete choice of \a k \ has to be taken with some care to optimize the con- 
vergence of the inverter. For example, taking \a k \ = 1 led to a bad convergence 
behaviour of the conjugate gradient algorithm. It is our experience that choos- 
ing \a k \ not much higher than \a r \, r being the index of the largest eigenvalue 
projected out, see eq. (|22]) , leads to a normal behaviour of the conjugate gradient 
algorithm. 

As an example of the effect of the improvement using M, eq. (|2"1~D, we show in 
Fig. ^| the behaviour of the pion propagator at zero distance T 7r (0) at (3 = 5.85 on 
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Figure 3: The pion propagator at zero distance as a function of N. The black 
squares are the results for the original 5D domain wall operator. The open sym- 
bols correspond to the improved operator with three and ten eigenvalues projected 
out. The data are obtained on an 8 3 • 16 lattice at f3 — 5.85 and am = 0.1. 
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an 8 3 ■ 16 lattice and for a quark mass, am = 0.1. A similar behaviour is obtained 
for the pion propagator at other distances. 

Already the projection of only three low-lying eigenvalues is sufficient to accelerate 
the convergence substantially: similar approximations to the N — > oo limit are 
obtained for N ~ 150 in the unimproved case and N ~ 30 in the improved one. It 
would, of course, be interesting to see the effect also on other physical quantities. 

Conclusions 



In this note we presented numerical evidence that in practical simulations domain 
wall fermions need an unacceptably large number of slices in the extra dimension 
to ensure that the chiral symmetry-breaking terms are suppressed. The reason is 
that very small eigenvalues of a 4D operator appear, which are directly related 
to the convergence rate of the 5D domain-wall operator. These results confirm 



the findings in |L9 



As in the case of Neuberger's operator, it is possible, however, to separate a 
number of eigenvalues of the 4D operator and treat them exactly or shift them 
to larger harmless values. We tested two different implementations of this idea. 
The first has already appeared in |jl9| , the second, which is described in detail 
above, is new. We found that numerically the second implementation works much 
better. 

It is our overall impression, however, that there is no particular advantage in 
using domain-wall fermions instead of Neuberger's operator. Theoretical consid- 
erations demonstrate that both approaches to realize a chiral symmetry on the 
lattice are equivalent, to the extent that they satisfy the Ginsparg- Wilson rela- 
tion. However, in practical applications, it is our present experience that it is 
easier to control chiral symmetry violations with Neuberger's operator. 
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